# some useful imports and settings
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from ipywidgets import interactive
import scipy.fft
import scipy.signal
import cv2
%config InlineBackend.figure_format = 'retina' # better looking figures on high-resolution screens
# automatically reload modules when running, otherwise jupyter does not notice if our functions have changed
%load_ext autoreload
%autoreload 2
# not used in this lab:
#import netCDF4 as nc
#import datetime
#import scipy.optimize
#import time
#import multiprocessing as mp
a. Read in the landsat data from 10/16/2023 over Albuquerque, and plot the bands separately.
The bands represent the following information:
b. Create an "RGB" composite to view the image in true color. Set the bounds of your plot to view just the area around Albuquerque and the Sandias.
c. Then try computing NDVI using the Red and Near-infrared bands and plot that image.
blue = cv2.imread('LC09_L2SP_033036_20231016_20231017_02_T1_SR_B2.TIF', cv2.IMREAD_UNCHANGED)
green = cv2.imread('LC09_L2SP_033036_20231016_20231017_02_T1_SR_B3.TIF')
red = cv2.imread('LC09_L2SP_033036_20231016_20231017_02_T1_SR_B4.TIF')
near_infrared = cv2.imread('LC09_L2SP_033036_20231016_20231017_02_T1_SR_B5.TIF')
# It seems that my version of OpenCV doesn't like reading TIF files by default and every solution I've found online
# doesn't seem to fix anything. My code, as well as the code in the notes file when I run it again, just gives
# a completely black image. My best guess is that TIF files are using a different system to record the brightness
# of each channel that OpenCV doesn't like to scale properly. I have written my code as though it does actually work,
# but of course it doesn't. I'm also guessing that you used an older version of OpenCV to get your code to read the exact
# same file properly. It's also possible that something went wrong in the download when these files were compressed!
fig, axes = plt.subplots(2,2,figsize=(10,10))
axes[0,0].imshow(blue)
axes[0,0].set_title('Blue')
axes[1,0].imshow(green)
axes[1,0].set_title('Green')
axes[0,1].imshow(red)
axes[0,1].set_title('Red')
axes[1,1].imshow(near_infrared)
axes[1,1].set_title('Near Infrared')
combined = cv2.merge((red, green, blue))
plt.imshow(combined)
plt.show()
ndvi = (near_infrared - red) / (near_infrared + red)
plt.imshow(ndvi)
plt.show()
"\nfig, axes = plt.subplots(2,2,figsize=(10,10))\naxes[0,0].imshow(blue)\naxes[0,0].set_title('Blue')\naxes[1,0].imshow(green)\naxes[1,0].set_title('Green')\naxes[0,1].imshow(red)\naxes[0,1].set_title('Red')\naxes[1,1].imshow(near_infrared)\naxes[1,1].set_title('Near Infrared')\n"
Apply the following filters to the grayscale image of my cat Pumpkin, and plot them all together with the original
For the last one, use the function 'cv2.canny', which is a common edge detection filter. You'll have to look up how to use it!
cat = cv2.imread('pumpkin.jpg')
cat_blur = cv2.blur(cat,(9,9))
cat_median = cv2.medianBlur(cat,9)
cat_gray = cv2.cvtColor(cat,cv2.COLOR_BGR2GRAY)
cat_gray_blur = cv2.GaussianBlur(cat_gray, (3, 3), 0)
cat_edge = cv2.Canny(image=cat_gray_blur, threshold1=100, threshold2=200)
fig,axs = plt.subplots(2,2,figsize=(20,20))
axs[0,0].imshow(cat)
axs[0,1].imshow(cat_blur)
axs[1,0].imshow(cat_median)
axs[1,1].imshow(cat_edge)
[155 187 138]
<matplotlib.image.AxesImage at 0x27bbb64f490>
The information output by the connected components identifier contains other useful information besides just the shape areas. For example, we can get the object centroids within our for loop, with
(cx, cy) = centroids[i]
and bounding-box information with
x = stats[i, cv2.CC_STAT_LEFT]
y = stats[i, cv2.CC_STAT_TOP]
w = stats[i, cv2.CC_STAT_WIDTH]
h = stats[i, cv2.CC_STAT_HEIGHT]
Copy the code above and modify it to plot the original grayscale image, but with a red text label on the figure at the centroid location of each iceberg displaying that iceberg's area if it's larger than your chosen threshold, and draw a bounding box around the largest one.
# load the (grayscale) image and make a plot
img = cv2.imread('iceye_icebergs.png')
img_gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
fig,axs = plt.subplots(1,2,figsize=(20,10))
axs[0].imshow(img_gray,cmap='gray')
# threshold the image at a grayscale value of 127 (halfway through the range).
# we also have to specify the "threshold type", which here is cv2.THRESH_BINARY
# finally, notice that there are 2 outputs of this function - the threshold value and the image.
# see a nice tutorial here: https://opencv24-python-tutorials.readthedocs.io/en/latest/py_tutorials/py_imgproc/py_thresholding/py_thresholding.html
thres_value, img_binary = cv2.threshold(img_gray,127,255,cv2.THRESH_BINARY)
axs[1].imshow(img_binary,cmap='gray')
plt.show()
# now apply the segmentation
# based (in part) on a long tutorial here:
# https://pyimagesearch.com/2021/02/22/opencv-connected-component-labeling-and-analysis/
# this value defines what is "connected" - it is either 4 (edges only) or 8 (diagonals OK)
connectivity = 4
# get the connected components:
numLabels, labels, stats, centroids = cv2.connectedComponentsWithStats(img_binary, connectivity, cv2.CV_32S)
# how many did we find?
print('identified', numLabels, 'icebergs')
# loop over the components, and if their areas are above the required size, add them to the map.
minArea = 1000
mask=np.zeros(np.shape(img_binary))
count=0
cx_list = []
cy_list = []
area_label_list = []
for i in range(numLabels):
area = stats[i,cv2.CC_STAT_AREA]
(cx, cy) = centroids[i]
# note, we skip label 1 because this is the largest - usually the background of the image!
if i>0 and area > minArea:
#componentMask = (labels == i).astype("uint8") * 255
#mask = cv2.bitwise_or(mask, componentMask)
count += 1
object_mask = (labels == i)
mask[object_mask] = area
cx_list.append(cx)
cy_list.append(cy)
area_label_list.append(area)
if i>0 and area >= 10000:
x = stats[i, cv2.CC_STAT_LEFT]
y = stats[i, cv2.CC_STAT_TOP]
w = stats[i, cv2.CC_STAT_WIDTH]
h = stats[i, cv2.CC_STAT_HEIGHT]
# make any areas that are still zero into nan, for a white background:
mask[mask==0] = np.nan
print(count, 'are larger than the min size of',minArea)
fig, ax = plt.subplots(1,1,figsize=(20,10))
ax.imshow(mask)
#fig.colorbar(cm.ScalarMappable(norm=norm, cmap=cmap), ax=ax)
rect = patches.Rectangle((x, y), w, h, linewidth=1, edgecolor='r', facecolor='none')
ax.add_patch(rect)
for i in range(len(cx_list)):
ax.text(cx_list[i],cy_list[i],area_label_list[i],
color="red",fontsize=15)
identified 2176 icebergs 7 are larger than the min size of 1000
Find an image you think could be easily thresholded and separated into connected components, then apply the above method to it. Be creative - for example, maybe you could use this method to identify and count balloons in the sky, or olivine crystals from a thin section, or even get the relative abundances of different minerals, if you can identify them separately using the different color bands.
Or maybe you can count penguins!
img = cv2.imread('balloons.jpg')
img_gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
fig,axs = plt.subplots(1,2,figsize=(20,10))
axs[0].imshow(img_gray,cmap='gray')
thres_value, img_binary = cv2.threshold(img_gray,135,255,cv2.THRESH_BINARY)
axs[1].imshow(img_binary,cmap='gray')
plt.show()
connectivity = 4
# get the connected components:
numLabels, labels, stats, centroids = cv2.connectedComponentsWithStats(img_binary, connectivity, cv2.CV_32S)
# how many did we find?
print('identified', numLabels, 'icebergs')
# loop over the components, and if their areas are above the required size, add them to the map.
minArea = 1000
mask=np.zeros(np.shape(img_binary))
count=0
cx_list = []
cy_list = []
area_label_list = []
for i in range(numLabels):
area = stats[i,cv2.CC_STAT_AREA]
(cx, cy) = centroids[i]
# note, we skip label 1 because this is the largest - usually the background of the image!
if i>0 and area > minArea:
#componentMask = (labels == i).astype("uint8") * 255
#mask = cv2.bitwise_or(mask, componentMask)
count += 1
object_mask = (labels == i)
mask[object_mask] = area
# make any areas that are still zero into nan, for a white background:
mask[mask==0] = np.nan
fig, ax = plt.subplots(1,1,figsize=(20,10))
ax.imshow(mask)
#fig.colorbar(cm.ScalarMappable(norm=norm, cmap=cmap), ax=ax)
identified 584 icebergs
<matplotlib.image.AxesImage at 0x2648f3d0490>